%--------------------------------------------------------------------------
% Script to plot likelihood densities
%--------------------------------------------------------------------------

clear; 
clc;
close all

set(0,'defaultTextInterpreter','latex');

LHPath = [pwd, '/Results/'];
NTstr  = 'N10000_T25';

% Load LogLH values
LogLH_PF1 = csvread( [LHPath, 'PFwithX_K8_M50/', NTstr, '_PFLogLH_M50.csv'],1,0);
LogLH_PF2 = csvread( [LHPath, 'PFwithX_K8_M500/', NTstr, '_PFLogLH_M500.csv'],1,0);
LogLH_PF3 = csvread( [LHPath, 'PFwithX_K8_M5000/', NTstr, '_PFLogLH_M5000.csv'],1,0);
LogLH_KF  = csvread( [LHPath, 'KFnoX_K8/', NTstr, '_lh_filtered.csv'],1,0);

% Rows t=1...T contain likelihood increments
% Row  t=T+1 contains sum of increments
[Tp1 Nsim] = size(LogLH_PF1);

% Load Nobs values
N_all = csvread( [LHPath, 'PFwithX_K8/', NTstr, '_N_all.csv'],1,0);
N_all = N_all + 2; %add aggregate observations
N_all(1,:) = zeros(size(N_all(1,:)));
N_total = sum(N_all);

% Select last row and standardize by the number of observations
LogLH_PF1 = (LogLH_PF1(Tp1,:)/N_total)';
LogLH_PF2 = (LogLH_PF2(Tp1,:)/N_total)';
LogLH_PF3 = (LogLH_PF3(Tp1,:)/N_total)';
LogLH_KF  = (LogLH_KF(Tp1,:)/N_total)';

% Compute mean of PF3 as benchmark
LogLH_true = log(mean(exp(LogLH_PF3)))

% Express in terms of deviations from KF and standardize by number of
% observations
LogLH_PF1 = LogLH_PF1 - LogLH_true;
LogLH_PF2 = LogLH_PF2 - LogLH_true;
LogLH_PF3 = LogLH_PF3 - LogLH_true;
LogLH_KF  = LogLH_KF  - LogLH_true;

%%
[LHdens_PF1_Y,LHdens_PF1_X] = ksdensity(LogLH_PF1);
[LHdens_PF2_Y,LHdens_PF2_X] = ksdensity(LogLH_PF2);
[LHdens_PF3_Y,LHdens_PF3_X] = ksdensity(LogLH_PF3);

ymin = 0;
ymax = 2E3;
xmin = -0.007;
xmax = 0.002;

figure(1);clf;
set(figure(1),'PaperType','usletter','PaperOrientation','Landscape','PaperPosition',[0.1 0.1 11 8.5]);
plot(LHdens_PF3_X,LHdens_PF3_Y,'Color','b','LineStyle','-','LineWidth',4)
hold on
plot(ones(size(LHdens_PF3_X))*LogLH_KF, linspace(ymin,ymax,length(LHdens_PF3_X)), ...
    'Color', 'k', 'LineStyle', '-', 'LineWidth',2)

set(gca,'FontSize',30)
% legend('M=50','M=500','M=5,000','Location','northeast')
axis([xmin xmax ymin ymax]) % for N=1,250
saveas(figure(1), [LHPath, NTstr, '_PFwithXLogLHdens.pdf'] );
%close all


